﻿// mod09ga_prepare.cpp : 此文件包含 "main" 函数。程序执行将在此处开始并结束。
//使用C盘 OSGeo4W的GDAL库快速开发GDAL程序。
//云遥服务器上 探测中心雷达拼图v3.0产品 处理tiff格式程序。 2025/4/16
// Z_RADA_C_BABJ_20250416010700_P_DOR_ACHN_CAP_20250416_010000_23.bin
// 将24层数据转为 500m 间隔的Geotiff格式数据。8km以上的通过插值得到500m间隔层。
// 结果总计32层 500m-16000m
// 层数对应关系
//           原始  结果
// 500-8000  0-15  0-15
// 8500      --    16
// 9000      16    17
// ...
// 15500     --    30
// 16000     23    31
//
// 32层数据量太大了，采用1000米间隔高度层，共15层
//           原始  结果
// 1000      1     0
// 2000      3,5,7,9,11,13     1,2,3,4,5,6
// 8000      15    7
// 9000      16    8
// 16000     23    15

#include <algorithm>
#include <iostream>
#include <string>
#include <vector>
#include <fstream>
using namespace std;
#include "gdal_priv.h"
#include "bzlib.h"
#include "wstringutils.h"
#define UNICODE
#include "tinydir.h"

struct Header {
    char label[4];
    char version[4];
    int file_bytes = 0;
    short mos_id = 0;
    short coordinate = 0;
    char varname[8];
    char desc[64];
    int blockPos = 0;
    int blockLen = 0;
    int timeZone = 0;
    short yr = 0;
    short mon = 0;
    short day = 0;
    short hr = 0;
    short minu = 0;
    short sec = 0;
    int obsSec = 0;
    short obsDates = 0;
    short genDates = 0;
    int genSec = 0;
    int latmin = 0;
    int lonmin = 0;
    int latmax = 0;
    int lonmax = 0;
    int cx = 0;
    int cy = 0;
    int xsize = 0;
    int ysize = 0;
    int dx = 0;
    int dy = 0;
    short height = 0;
    short compress = 0;
    int numradars = 0;
    int unzipbytes = 0;
    short scale = 0;
    short unused = 0;
    char rgnid[8];
    char units[8];
    char reserved[60];

};

/// 结果会按照文件名字符串从小到大排序。
/// 结果只包括文件名，没有路径信息。
vector<string> scanFiles(string indir, string dateTimeStr) {
    tinydir_dir dir;
    wstring wdirstr = wstring(indir.begin(), indir.end());
    wstring wdtStr = wstring(dateTimeStr.begin(), dateTimeStr.end());
    tinydir_open(&dir, wdirstr.c_str());
    vector<string> res;
    while (dir.has_next)
    {
        tinydir_file file;
        tinydir_readfile(&dir, &file);
        wstring fname1(  file.name);
        if (fname1.find(wdtStr) != string::npos) {
            res.push_back( string(fname1.begin(),fname1.end() ));
        }
        tinydir_next(&dir);
    }
    tinydir_close(&dir);
    std::sort(res.begin(), res.end());
    return res;
}
//筛选出第一个处理时间
vector<string> filterFirstProcess(vector<string> fnameArr) {
    // filename should be like Z_RADA_C_BABJ_20250415013103_P_DOR_ACHN_CAP_20250415_012400_07.bin
    vector<string> res;
    if (fnameArr.size() > 0) {
        string filename0 = fnameArr[0];
        if (filename0.size() == 66 || filename0.size()==64 ) {//64 for cref
            string processDtStr = filename0.substr(14,14);
            for (int i = 0; i < fnameArr.size(); ++i) {
                if (fnameArr[i].find(processDtStr) != string::npos) {
                    res.push_back(fnameArr[i]);
                }
            }
        }
    }
    return res;
}


/// dateTimeStr should be 20250416_010000 format.
/// return 0 is OK.
///  the first part of datetime maybe has two different processing time. use the first process time.
int process32Layers(string inDayDir, string dateTimeStr, string outfile) {

    vector<string> filenameArr = scanFiles(inDayDir, dateTimeStr);
    filenameArr = filterFirstProcess(filenameArr);
    if (filenameArr.size() != 24) {
        std::cout << "file count not 24" << endl;
        return 12;
    }
    CPLStringList opt;
    opt.SetNameValue("COMPRESS", "LZW");
    const int outxsize = 6200;
    const int outysize = 4200;
    GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
    GDALDataset* dsout = driver->Create(outfile.c_str(),
        outxsize, outysize, 32, GDT_Int16, opt.List());
    if (dsout == 0) {
        std::cout << "failed to create output " << outfile << endl;
        return 13;
    }
    vector<unsigned char> lastDatabuffer(outxsize * outysize * 2);
    int bandIndex0Based = 0;
    for (int i = 0; i < filenameArr.size() ; ++i) {

        string filenameX = inDayDir+"/"+filenameArr[i] ;
        //sprintf(filename, "E:/TempAndClean/Z_RADA_C_BABJ_20250416010700_P_DOR_ACHN_CAP_20250416_010000_%02d.bin", i);
        FILE* pf = fopen(filenameX.c_str(), "rb");
        if (pf == 0) {
            std::cout << "failed to open " << filenameX << endl;
            return 14;
        }

        Header header1;
        fread(&header1, sizeof(Header), 1, pf);
        int     bzerror = 0;
        if (header1.xsize != outxsize) {
            std::cout << "wrong header1.xsize " << header1.xsize << endl;
            return 15;
        }
        if (header1.ysize != outysize) {
            std::cout << "wrong header1.ysize " << header1.xsize << endl;
            return 16;
        }
        //the pf is already move to 256 , so no need fseak.
        //fseek(pf, 256, SEEK_SET);
        BZFILE* bf = BZ2_bzReadOpen(&bzerror, pf, 0, 0, NULL, 0);
        if (bzerror != BZ_OK) {
            BZ2_bzReadClose(&bzerror, bf); bf = 0;
            fclose(pf); pf = 0;
            std::cout << "BZ2_bzReadOpen is failed." << endl;
            return 17;
        }
        vector<unsigned char> databuffer(header1.xsize * header1.ysize * 2);//data is short.
        int nread = BZ2_bzRead(&bzerror, bf, databuffer.data(), databuffer.size());
        if (bzerror != BZ_OK && bzerror != BZ_STREAM_END) {
            //BZ_STREAM_END is OK, not error.
            std::cout << "bad  " << bzerror << endl;
        }
        if (bf != 0) BZ2_bzReadClose(&bzerror, bf); bf = 0;
        if (pf != 0) fclose(pf); pf = 0;

        int asize = outxsize * outysize;
        short* sdata = (short*)databuffer.data();
        for (int it = 0; it < asize; ++it) {
            if (sdata[it] < 0) {
                sdata[it] = 0;//use 0 to fill all negative values.
            }
        }
        std::cout << i << " height " << header1.height << " scale " << header1.scale << endl;
        if (i <= 15) {
            dsout->GetRasterBand(bandIndex0Based + 1)->RasterIO(GF_Write, 0, 0, outxsize, outysize, databuffer.data(), outxsize, outysize, GDT_Int16, 0, 0, 0);
            bandIndex0Based++;
        }
        else {
            //interpol between heights.
            short* dataLow = (short*)lastDatabuffer.data();
            short* dataHigh = (short*)databuffer.data();
            for (int it = 0; it < asize; ++it) {
                if (dataLow[it] > 0 && dataHigh[it] > 0) {
                    dataLow[it] = (dataLow[it] + dataHigh[it]) / 2.0;
                }
                else {
                    dataLow[it] = 0;//the value lower than 0 are all fill with 0.
                }
            }
            dsout->GetRasterBand(bandIndex0Based + 1)->RasterIO(GF_Write, 0, 0,
                outxsize, outysize, lastDatabuffer.data(), outxsize, outysize, GDT_Int16, 0, 0, 0);
            bandIndex0Based++;
            dsout->GetRasterBand(bandIndex0Based + 1)->RasterIO(GF_Write, 0, 0,
                outxsize, outysize, databuffer.data(), outxsize, outysize, GDT_Int16, 0, 0, 0);
            bandIndex0Based++;
        }
        std::cout << "bandIndex0Based " << bandIndex0Based << endl;
        if (i >= 15) {
            lastDatabuffer = databuffer;
        }
    }
    double trans[] = { 73.0, 0.01, 0.00 , 54.2, 0.00, -0.01 };
    dsout->SetGeoTransform(trans);     // trans[6]
    dsout->SetProjection("GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],"
        "AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,"
        "AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]");     //使用这个不需要注册环境变量，也不需要引用ogr的头文件
    GDALClose(dsout);
    return 0;
}

// 1000      1     0
// 2000      3,5,7,9,11,13     1,2,3,4,5,6
// 8000      15    7
// 9000      16    8
// 16000     23    15
int process16Layers(string inDayDir, string dateTimeStr, string outfile) {

    vector<string> filenameArr = scanFiles(inDayDir, dateTimeStr);
    filenameArr = filterFirstProcess(filenameArr);
    if (filenameArr.size() != 24) {
        std::cout << "file count not 24" << endl;
        return 12;
    }
    CPLStringList opt;
    opt.SetNameValue("COMPRESS", "LZW");
    const int outxsize = 6200;
    const int outysize = 4200;
    GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
    GDALDataset* dsout = driver->Create(outfile.c_str(),
        outxsize, outysize, 16, GDT_Int16, opt.List());
    vector<int> usedLayers = { 1,3,5,7,9,11,13,15,16,17,18,19,20,21,22,23 };
    if (dsout == 0) {
        std::cout << "failed to create output " << outfile << endl;
        return 13;
    }
    int bandIndex0Based = 0;
    for (int i = 0; i < filenameArr.size(); ++i) {
        bool isUse = false;
        for (int j = 0; j < usedLayers.size(); ++j) {
            if (usedLayers[j] == i) {
                isUse = true;
                break;
            }
        }
        if (isUse == false)continue;//skip
        string filenameX = inDayDir + "/" + filenameArr[i];
        //sprintf(filename, "E:/TempAndClean/Z_RADA_C_BABJ_20250416010700_P_DOR_ACHN_CAP_20250416_010000_%02d.bin", i);
        FILE* pf = fopen(filenameX.c_str(), "rb");
        if (pf == 0) {
            std::cout << "failed to open " << filenameX << endl;
            return 14;
        }

        Header header1;
        fread(&header1, sizeof(Header), 1, pf);
        int     bzerror = 0;
        if (header1.xsize != outxsize) {
            std::cout << "wrong header1.xsize " << header1.xsize << endl;
            return 15;
        }
        if (header1.ysize != outysize) {
            std::cout << "wrong header1.ysize " << header1.xsize << endl;
            return 16;
        }
        //the pf is already move to 256 , so no need fseak.
        //fseek(pf, 256, SEEK_SET);
        BZFILE* bf = BZ2_bzReadOpen(&bzerror, pf, 0, 0, NULL, 0);
        if (bzerror != BZ_OK) {
            BZ2_bzReadClose(&bzerror, bf); bf = 0;
            fclose(pf); pf = 0;
            std::cout << "BZ2_bzReadOpen is failed." << endl;
            return 17;
        }
        vector<unsigned char> databuffer(header1.xsize * header1.ysize * 2);//data is short.
        int nread = BZ2_bzRead(&bzerror, bf, databuffer.data(), databuffer.size());
        if (bzerror != BZ_OK && bzerror != BZ_STREAM_END) {
            //BZ_STREAM_END is OK, not error.
            std::cout << "bad  " << bzerror << endl;
        }
        if (bf != 0) BZ2_bzReadClose(&bzerror, bf); bf = 0;
        if (pf != 0) fclose(pf); pf = 0;

        int asize = outxsize * outysize;
        short* sdata = (short*)databuffer.data();
        for (int it = 0; it < asize; ++it) {
            if (sdata[it] < 0) {
                sdata[it] = 0;//use 0 to fill all negative values.
            }
        }
        std::cout << i << " height " << header1.height << " scale " << header1.scale << endl;
        dsout->GetRasterBand(bandIndex0Based + 1)->RasterIO(GF_Write, 0, 0, outxsize, outysize,
            databuffer.data(), outxsize, outysize, GDT_Int16, 0, 0, 0);
        bandIndex0Based++;
    }
    double trans[] = { 73.0, 0.01, 0.00 , 54.2, 0.00, -0.01 };
    dsout->SetGeoTransform(trans);     // trans[6]
    dsout->SetProjection("GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],"
        "AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,"
        "AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]");     //使用这个不需要注册环境变量，也不需要引用ogr的头文件
    GDALClose(dsout);
    return 0;
}

//one layer
int processCref(string inDayDir, string dateTimeStr, string outfile) {

    vector<string> filenameArr = scanFiles(inDayDir, dateTimeStr);
    filenameArr = filterFirstProcess(filenameArr);
    if (filenameArr.size() == 0) {
        std::cout << "file count == 0" << endl;
        return 12;
    }
    CPLStringList opt;
    opt.SetNameValue("COMPRESS", "LZW");
    const int outxsize = 6200;
    const int outysize = 4200;
    GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
    GDALDataset* dsout = driver->Create(outfile.c_str(),
        outxsize, outysize, 1, GDT_Int16, opt.List());
    if (dsout == 0) {
        std::cout << "failed to create output " << outfile << endl;
        return 13;
    }

    string filenameX = inDayDir + "/" + filenameArr[0];
    //sprintf(filename, "E:/TempAndClean/Z_RADA_C_BABJ_20250416010700_P_DOR_ACHN_CAP_20250416_010000_%02d.bin", i);
    FILE* pf = fopen(filenameX.c_str(), "rb");
    if (pf == 0) {
        std::cout << "failed to open " << filenameX << endl;
        return 14;
    }

    Header header1;
    fread(&header1, sizeof(Header), 1, pf);
    int     bzerror = 0;
    if (header1.xsize != outxsize) {
        std::cout << "wrong header1.xsize " << header1.xsize << endl;
        return 15;
    }
    if (header1.ysize != outysize) {
        std::cout << "wrong header1.ysize " << header1.xsize << endl;
        return 16;
    }
    //the pf is already move to 256 , so no need fseak.
    //fseek(pf, 256, SEEK_SET);
    BZFILE* bf = BZ2_bzReadOpen(&bzerror, pf, 0, 0, NULL, 0);
    if (bzerror != BZ_OK) {
        BZ2_bzReadClose(&bzerror, bf); bf = 0;
        fclose(pf); pf = 0;
        std::cout << "BZ2_bzReadOpen is failed." << endl;
        return 17;
    }
    vector<unsigned char> databuffer(header1.xsize * header1.ysize * 2);//data is short.
    int nread = BZ2_bzRead(&bzerror, bf, databuffer.data(), databuffer.size());
    if (bzerror != BZ_OK && bzerror != BZ_STREAM_END) {
        //BZ_STREAM_END is OK, not error.
        std::cout << "bad  " << bzerror << endl;
    }
    if (bf != 0) BZ2_bzReadClose(&bzerror, bf); bf = 0;
    if (pf != 0) fclose(pf); pf = 0;

    int asize = outxsize * outysize;
    short* sdata = (short*)databuffer.data();
    for (int it = 0; it < asize; ++it) {
        if (sdata[it] < 0) {
            sdata[it] = 0;//use 0 to fill all negative values.
        }
    }
    dsout->GetRasterBand(1)->RasterIO(GF_Write, 0, 0, outxsize, outysize,
        databuffer.data(), outxsize, outysize, GDT_Int16, 0, 0, 0);

    double trans[] = { 73.0, 0.01, 0.00 , 54.2, 0.00, -0.01 };
    dsout->SetGeoTransform(trans);     // trans[6]
    dsout->SetProjection("GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],"
        "AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,"
        "AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]");     //使用这个不需要注册环境变量，也不需要引用ogr的头文件
    GDALClose(dsout);
    return 0;
}

struct DataHeader {
    byte pcode;
    int scale;
    int offset;
    char binlen;
    char flags;
    char ttt;
    //short rowres;//meter
    //short colres;//meter
    //int hgrids;
    //int vgrids;
    //float maxvalue;
    //int hcmax;//最大值横纵坐标
    //int vcmax;
    //float minvalue;
    //int hcmin;
    //int vcmin;
    //char reserve[17];

};


const float D2R = 3.1415926 / 180.0;

void sd2uv(float spd, float deg, float& u, float& v) {
    u = -9999.0;
    v = -9999.0;
    if (spd > -9998.f && deg > -9998.0) {
        u = spd * cosf(D2R * deg);
        v = spd * sinf(D2R * deg);
    }
}


/// <summary>
/// 92 是水平风速 ， 93是0-360风向。
/// 36k reso 150x150 size 94x94
/// 96k reso 90x90   size 156x156
/// 629k reso 50x50  size 401x401
/// </summary>
/// <param name="filename92"> 92 是水平风速 </param>
/// <param name="outfile"></param>
/// <returns>返回0成功，其他失败</returns>
int processRadarWind(string filename92, string filename93,  string outpngfilename, string outjsonfilename) {

    //sprintf(filename, "E:/TempAndClean/Z_RADA_C_BABJ_20250416010700_P_DOR_ACHN_CAP_20250416_010000_%02d.bin", i);
    FILE* pf = fopen(filename92.c_str(), "rb");
    if (pf == 0) {
        std::cout << "failed to open " << filename92 << endl;
        return 14;
    }

    FILE* pf93 = fopen(filename93.c_str(), "rb");

    char header0[32];
    char header1[128];
    char header2[256];
    char header3[128];
    char header4[64];
    char header5[64];
    fread(header0 , 32 , 1 , pf );
    fread(header1, 128, 1, pf);
    fread(header2, 256, 1, pf);
    fread(header3, 128, 1, pf);
    fread(header4, 64, 1, pf);
    fread(header5, 64, 1, pf); 
    int scantype92 = (int)header2[160];
    int cutnumber92 = (int)header2[173];
    cout << "92 header2 scan type " << scantype92 << endl;
    cout << "92 header2 cut number " << cutnumber92 << endl;
    if (scantype92 != 2 && cutnumber92 != 1) {
        cout << "092 file scan type or cut number not supported. only support scan type=2 and cut number=1." << endl;
        return 15;
    }
    
    fread(header0, 32, 1, pf93);
    fread(header1, 128, 1, pf93);
    fread(header2, 256, 1, pf93);
    fread(header3, 128, 1, pf93);
    fread(header4, 64, 1, pf93);
    fread(header5, 64, 1, pf93);
    int scantype93 = (int)header2[160];
    int cutnumber93 = (int)header2[173];
    cout << "93 header2 scan type " << scantype93 << endl;
    cout << "93 header2 cut number " << cutnumber93 << endl;
    if (scantype93 != 2 && cutnumber93 != 1) {
        cout << "093 file scan type or cut number not supported. only support scan type=2 and cut number=1." << endl;
        return 15;
    }

    char pcode;
    fread(&pcode, 1, 1, pf);
    cout << "92 pcode " << (int)pcode << endl;
    fread(&pcode, 1, 1, pf93);
    cout << "93 pcode " << (int)pcode << endl;
    int scale=0;
    int offset = 0;
    char binlen = 0;
    char flags = 0;
    fread(&scale, 4, 1, pf);
    fread(&offset, 4, 1, pf);
    fread(&binlen, 1, 1, pf);
    fread(&flags, 1, 1, pf);
    cout << "92 scale offset " << scale << " " << offset << endl;
    cout << "92 binLen flags " << (int)binlen << " " << (int)flags << endl;

    fread(&scale, 4, 1, pf93);
    fread(&offset, 4, 1, pf93);
    fread(&binlen, 1, 1, pf93);
    fread(&flags, 1, 1, pf93);

    cout << "93 scale offset " << scale << " " << offset << endl;
    cout << "93 binLen flags " << (int)binlen << " " << (int)flags << endl;

    short resox,resoy;
    fread(&resox, 2, 1, pf);
    fread(&resoy, 2, 1, pf);
    cout << "92 resox resoy " << resox << " " << resoy << endl;
    fread(&resox, 2, 1, pf93);
    fread(&resoy, 2, 1, pf93);
    cout << "93 resox resoy " << resox << " " << resoy << endl;

    int xsize = 0;
    int ysize = 0;
    fread(&xsize, 4, 1, pf);
    fread(&ysize, 4, 1, pf);
    cout << "92 xsize ysize " << xsize << " " << ysize << endl;
    fread(&xsize, 4, 1, pf93);
    fread(&ysize, 4, 1, pf93);
    cout << "93 xsize ysize " << xsize << " " << ysize << endl;

    const int OTHERLEN = 4 + 4 + 4 + 4 + 4 + 4 + 17;
    char otherHeaderData[OTHERLEN];
    fread(otherHeaderData, OTHERLEN, 1, pf);
    fread(otherHeaderData, OTHERLEN, 1, pf93);
    //92风速
    //93风向

    vector<float> dataSpd0(xsize * ysize);
    vector<float> dataDeg0(xsize * ysize);
    int it = 0;
    for (int iy = 0; iy < ysize; ++iy) {
        for (int ix = 0; ix < xsize; ++ix) {
            float vvv;
            fread(&vvv, 4, 1, pf);
            dataSpd0[it] = vvv;
            fread(&vvv, 4, 1, pf93);
            dataDeg0[it] = vvv;
            //if (vvv > -9998) {
            //    cout << vvv << endl;
            //}
            it++;
        }
    }
    fclose(pf93);
    fclose(pf);

    //无论输入分辨率和尺寸多少，输出统一使用分辨率50，尺寸401
    const int outxsize = 401;
    const int outysize = 401;
    const float outresox = 50;
    const float outresoy = 50;
    
    //重采样到401x401
    vector<float> dataSpd(outxsize * outysize);
    vector<float> dataDeg(outxsize * outysize);
    float outx0 = outxsize / 2.f;
    float outy0 = outysize / 2.f;
    float inx0 = xsize / 2.f;
    float iny0 = ysize / 2.f;
    for (int iy = 0; iy < outysize; ++iy) {
        for (int ix = 0; ix < outxsize; ++ix) {
            float outcoorx = (ix - outx0) * outresox;
            float outcoory = (iy - outy0) * outresoy;
            int inx = outcoorx/resox + inx0 ;
            int iny = outcoory/resoy + iny0 ;
            if (inx >= 0 && inx < xsize && iny >= 0 && iny < ysize) {
                dataSpd[iy * outxsize + ix] = dataSpd0[iny * xsize + inx];
                dataDeg[iy * outxsize + ix] = dataDeg0[iny * xsize + inx];
            }
            else {
                dataSpd[iy * outxsize + ix] = -9999;
                dataDeg[iy * outxsize + ix] = -9999;
            }


        }
    }



    {
        vector<unsigned char> r(outxsize * outysize);
        vector<unsigned char> g(outxsize * outysize);
        vector<unsigned char> b(outxsize * outysize);
        vector<unsigned char> a(outxsize * outysize);



        /*
        0, 0, 255     2.5
48, 48, 255   5.0
101, 101, 255 10
179, 180, 255 15
230, 255, 230 20
179, 255, 179 25
101, 255, 101 30
0, 255, 0     35
        */
        int asize = outxsize * outysize;
        for (int it = 0; it < asize; ++it) {
            float spd = dataSpd[it];
            a[it] = 255;
            if (spd < -9998.f) {
                r[it] = 0;
                g[it] = 0;
                b[it] = 0;
                a[it] = 0;
            }
            else if (spd < 2.5) {
                r[it] = 0;
                g[it] = 0;
                b[it] = 255;
            }
            else if (spd < 5.0) {
                r[it] = 48;
                g[it] = 48;
                b[it] = 255;
            }
            else if (spd < 10.0) {
                r[it] = 101;
                g[it] = 101;
                b[it] = 255;
            }
            else if (spd < 15.0) {
                r[it] = 178;
                g[it] = 180;
                b[it] = 255;
            }
            else if (spd < 20.0) {
                r[it] = 230;
                g[it] = 255;
                b[it] = 230;
            }
            else if (spd < 25) {
                r[it] = 180;
                g[it] = 255;
                b[it] = 180;
            }
            else if (spd < 30) {
                r[it] = 101;
                g[it] = 255;
                b[it] = 101;
            }
            else {
                r[it] = 0;
                g[it] = 255;
                b[it] = 0;
            }
        }

        //
        GDALDriver* memDriver = GetGDALDriverManager()->GetDriverByName("MEM");
        GDALDataset* memds = memDriver->Create("", outxsize,
            outysize, 4, GDT_Byte, 0);
        memds->GetRasterBand(1)->RasterIO(GF_Write, 0, 0, outxsize, outysize, r.data(), outxsize, outysize, GDT_Byte, 0, 0, 0);
        memds->GetRasterBand(2)->RasterIO(GF_Write, 0, 0, outxsize, outysize, g.data(), outxsize, outysize, GDT_Byte, 0, 0, 0);
        memds->GetRasterBand(3)->RasterIO(GF_Write, 0, 0, outxsize, outysize, b.data(), outxsize, outysize, GDT_Byte, 0, 0, 0);
        memds->GetRasterBand(4)->RasterIO(GF_Write, 0, 0, outxsize, outysize, a.data(), outxsize, outysize, GDT_Byte, 0, 0, 0);
        memds->FlushCache();


        GDALDriver* d = GetGDALDriverManager()->GetDriverByName("PNG");

        GDALDataset* pngds = d->CreateCopy(outpngfilename.c_str(), memds, FALSE, 0, 0, 0);
        if (pngds == 0) {
            cout << "error, failed to create png output " << outpngfilename << endl;
            return 14;
        }

        GDALClose(memds);
        GDALClose(pngds);


    }

    {
        /*
        117.688350	-	117.888850  39.059450	-	39.259950
        */

        string part1 = "[{\"header\":{\"parameterCategory\":1,\"parameterNumber\":2,\"la1\":39.259950,\"la2\":39.059450,\"lo1\":117.688350,\"lo2\":117.888850,\"extent\":[117.688350,39.059450,117.888850,39.259950],\"nx\":401,\"ny\":401,\"dx\":1,\"dy\":1,\"min\":-20.0,\"max\":30.0,\"GRIB_COMMENT\":\"u-component of wind [m/s]\",\"GRIB_DISCIPLINE\":\"0(Meteorological)\",\"GRIB_ELEMENT\":\"UGRD\",\"GRIB_FORECAST_SECONDS\":\"0 sec\",\"GRIB_IDS\":\"CENTER=7(US-NCEP) SUBCENTER=0 MASTER_TABLE=2 LOCAL_TABLE=1 SIGNF_REF_TIME=1(Start_of_Forecast) REF_TIME=2020-06-20T00:00:00Z PROD_STATUS=0(Operational) TYPE=1(Forecast)\",\"GRIB_PDS_PDTN\":\"0\",\"GRIB_PDS_TEMPLATE_ASSEMBLED_VALUES\":\"2 2 2 0 81 0 0 1 0 103 0 10 255 0 0\",\"GRIB_PDS_TEMPLATE_NUMBERS\":\"2 2 2 0 81 0 0 0 1 0 0 0 0 103 0 0 0 0 10 255 0 0 0 0 0\",\"GRIB_REF_TIME\":\"1592611200 sec UTC\",\"GRIB_SHORT_NAME\":\"10-HTGL\",\"GRIB_UNIT\":\"[m/s]\",\"GRIB_VALID_TIME\":\"1592611200 sec UTC\"},\"data\":[";
        string part2 = "] }, { \"header\":{\"parameterCategory\":1,\"parameterNumber\":3,\"la1\":39.259950,\"la2\":39.059450,\"lo1\":117.688350,\"lo2\":117.88885,\"extent\":[117.688350,39.059450,117.88885,39.259950],\"nx\":401,\"ny\":401,\"dx\":1,\"dy\":1,\"min\":-20.0,\"max\":30.0,\"GRIB_COMMENT\":\"v-component of wind [m/s]\",\"GRIB_DISCIPLINE\":\"0(Meteorological)\",\"GRIB_ELEMENT\":\"VGRD\",\"GRIB_FORECAST_SECONDS\":\"0 sec\",\"GRIB_IDS\":\"CENTER=7(US-NCEP) SUBCENTER=0 MASTER_TABLE=2 LOCAL_TABLE=1 SIGNF_REF_TIME=1(Start_of_Forecast) REF_TIME=2020-06-20T00:00:00Z PROD_STATUS=0(Operational) TYPE=1(Forecast)\",\"GRIB_PDS_PDTN\":\"0\",\"GRIB_PDS_TEMPLATE_ASSEMBLED_VALUES\":\"2 3 2 0 81 0 0 1 0 103 0 10 255 0 0\",\"GRIB_PDS_TEMPLATE_NUMBERS\":\"2 3 2 0 81 0 0 0 1 0 0 0 0 103 0 0 0 0 10 255 0 0 0 0 0\",\"GRIB_REF_TIME\":\"1592611200 sec UTC\",\"GRIB_SHORT_NAME\":\"10-HTGL\",\"GRIB_UNIT\":\"[m/s]\",\"GRIB_VALID_TIME\":\"1592611200 sec UTC\"},\"data\":[";
        string part3 = "] }]";
        int asize = outxsize * outysize;
        char buff[32];

        vector<float> udata(asize);
        vector<float> vdata(asize);

        vector<unsigned char> utexdata(asize); //UV 风纹理数据 整数部分+150 
        vector<unsigned char> vtexdata(asize);
        vector<unsigned char> u2texdata(asize); //UV风小数部分*100
        vector<unsigned char> v2texdata(asize);

        float u = -9999.0;
        float v = -9999.0;
        for (int it = 0; it < asize; ++it) {
            sd2uv(dataSpd[it], dataDeg[it], u, v);
            if (u < -9998.0 || v < -9998.0) {
                u = 0.001; //无效值给一个很小值
                v = 0.001;

                utexdata[it] = 150;
                vtexdata[it] = 150;
                u2texdata[it] = 0;
                v2texdata[it] = 0;
            }
            else {
                utexdata[it] = u+150 ; // 将UV风数据从-100~+100 平移150，转换后Byte值为 50~250。使用 0 作为填充值。
                vtexdata[it] = v+150 ;
                u2texdata[it] = ( (u + 150) - (int)(u+150) ) * 100;
                v2texdata[it] = ((v + 150) - (int)(v + 150)) * 100;
            }
            udata[it] = u;
            vdata[it] = v;
        }
        
        // cesium Wind中会对无效值（不在min max范围内）的值取最小值

        sprintf(buff, "%8.4f", udata[0] );
        part1 += string(buff);
        sprintf(buff, "%8.4f", vdata[0]);
        part2 += string(buff);

        for (int it = 1; it < asize; ++it) {
            sprintf(buff, ",%8.4f", udata[it]);
            part1 += string(buff);
            sprintf(buff, ",%8.4f", vdata[it]);
            part2 += string(buff);

        }
        ofstream ofs(outjsonfilename.c_str());
        if (ofs.is_open() == false) {
            cout << "error, failed to create output json file " << outjsonfilename << endl;
            return 13;
        }
        ofs << part1 << part2 << part3;
        ofs.close();

        string outuvpngfilename = outpngfilename;
        outuvpngfilename = wStringUtils::replaceString(outuvpngfilename, ".png", ".UVTEX.png");

        GDALDriver* memDriver = GetGDALDriverManager()->GetDriverByName("MEM");
        GDALDataset* memds = memDriver->Create("", outxsize,
            outysize, 4, GDT_Byte, 0);
        memds->GetRasterBand(1)->RasterIO(GF_Write, 0, 0, outxsize, outysize, utexdata.data(), outxsize, outysize, GDT_Byte, 0, 0, 0);   //R
        memds->GetRasterBand(2)->RasterIO(GF_Write, 0, 0, outxsize, outysize, vtexdata.data(), outxsize, outysize, GDT_Byte, 0, 0, 0);   //G
        memds->GetRasterBand(3)->RasterIO(GF_Write, 0, 0, outxsize, outysize, u2texdata.data(), outxsize, outysize, GDT_Byte, 0, 0, 0);   //G
        memds->GetRasterBand(4)->RasterIO(GF_Write, 0, 0, outxsize, outysize, v2texdata.data(), outxsize, outysize, GDT_Byte, 0, 0, 0);//A
        memds->FlushCache();
        GDALDriver* d = GetGDALDriverManager()->GetDriverByName("PNG");
        GDALDataset* pngds = d->CreateCopy(outuvpngfilename.c_str(), memds, FALSE, 0, 0, 0);
        if (pngds == 0) {
            cout << "error, failed to create png output " << outpngfilename << endl;
            return 14;
        }
        cout << "write uv wind texture OK:" << outuvpngfilename << endl;
        GDALClose(memds);
        GDALClose(pngds);

    }

    return 0;
}



int main(int argc,char* argv[])
{
    std::cout << "v0.2 401x401 ok" << endl;
    std::cout << "v0.3 for 7km data" << endl;
    cout << "v0.4 fix value scale bug." << endl;
    cout << "v0.5 fix not reading 093 data block header bug." << endl;
    cout << "v0.6 output uv wind texture." << endl;
    std::cout << "usage:yy_park_radar /some/dir/...092.BIN /some/dir/...093.BIN /out/dir" << endl;
    if (argc != 4) {
        cout << "error, argc not 4." << endl;
        return 11;
    }

    GDALAllRegister();
    string file92 = argv[1];// "Z_RADR_I_H0327_20250310084653_P_CDWL_LR_092.BIN";
    string file93 = argv[2]; //"Z_RADR_I_H0327_20250310084653_P_CDWL_LR_093.BIN";
    string outdir = argv[3];
    size_t pos1 = file92.find("Z_RADR_I_");
    if (pos1 == string::npos) {
        cout << "error, input filename does not has 'Z_RADR_I_'." << endl;
        return 12;
    }

    string baseName = file92.substr(pos1);
    string outpng = outdir + "/" + wStringUtils::replaceString( baseName,"_092.BIN",".png") ;// "Z_RADR_I_H0327_20250310084653_P_CDWL_LR.png";
    string outjson = outdir + "/" + wStringUtils::replaceString(baseName, "_092.BIN", ".json");// "Z_RADR_I_H0327_20250310084653_P_CDWL_LR.json";
    int resultCode = processRadarWind(file92, file93,outpng,outjson);

    if( resultCode==0 ) std::cout << "done" << endl;
    return resultCode;
}

